Set up workspace, load data, and load required packages.
rm(list=ls(all=TRUE))
results <- read.csv("Data/Urchin_ClimateChange_Data.csv", header=T, na.strings="NA")
library("lme4") #linear mixed modeling
library("lmtest") #linear mixed modeling
library("lmerTest") #calculate p-values in models
library("car") #levenes tests
library("emmeans") #post-hoc tests
library("lsmeans") #post-hoc tests
library("plotrix") #calculate std.error
library("plyr") #splitting, applying, and combining data
library("dplyr")
library("cowplot") #grid plotting and arrange plots
library("effects") #plot effects of modeling
library("MuMIn")
library("car") #levenes tests, Anova, qqPlots,
library("tidyverse")
library("stats")
library("onewaytests") #allows for Welch test with unequal variances (spine lengths)
library("stats")
library("effsize")
library("FSA")
library("rcompanion")
library("broom") #turns model output into data frames
library("tidyverse")
1) Build and run model
2) Check for normality of residuals
3) Check for homogeniety of variance of residuals
4) Look at model summary
5) Run anova to check for significance of main effects
EnviroMean <-
results %>%
select("Day", "Temperature", "pH", "Temp","pHspec") %>%
filter(Day >= "1") %>%
drop_na(Day) %>%
group_by(Temperature, pH) %>%
summarise(meanT = mean(Temp, na.rm = T)); EnviroMean
## # A tibble: 4 x 3
## # Groups: Temperature [2]
## Temperature pH meanT
## <fct> <fct> <dbl>
## 1 ambient acidified 25.0
## 2 ambient ambient 25.0
## 3 heated acidified 27.1
## 4 heated ambient 27.0
Environmental <-
results %>%
select("Day", "Temp", "pH", "Temperature","pHspec","Header", "TankID") %>%
filter(Day >= "1", Day != "126") #use only days from start of experiment when desired conditions were reached
#temp between high and ambient temp
summary(aov(Temp ~ Temperature, data=Environmental)) #temperature was different between high and ambient treatment
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 1 782.3 782.3 307.6 <2e-16 ***
## Residuals 771 1960.6 2.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 355 observations deleted due to missingness
#pH between high and ambient co2
summary(aov(pHspec ~ pH, data=Environmental)) #pH was different between high and ambient co2 treatments
## Df Sum Sq Mean Sq F value Pr(>F)
## pH 1 14.012 14.012 5247 <2e-16 ***
## Residuals 627 1.674 0.003
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 499 observations deleted due to missingness
##HEADERS
#temp between headers in high
aov1<- Environmental %>%
filter(Temperature=="heated")
summary(aov(Temp ~ as.factor(Header), data=aov1)) #not different between headers
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Header) 3 2.0 0.674 0.277 0.842
## Residuals 381 925.7 2.430
## 179 observations deleted due to missingness
#temp between headers in ambient temp
aov2<- Environmental %>%
filter(Temperature=="ambient")
summary(aov(Temp ~ as.factor(Header), data=aov2)) #not different between headers
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Header) 3 0.1 0.0318 0.012 0.998
## Residuals 384 1032.8 2.6897
## 176 observations deleted due to missingness
#pH between headers in acidified
aov3<- Environmental %>%
filter(pH=="acidified")
summary(aov(pHspec ~ as.factor(Header), data=aov3)) #not different between headers
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Header) 3 0.0243 0.008111 1.785 0.15
## Residuals 316 1.4360 0.004544
## 248 observations deleted due to missingness
#pH between headers in ambient pH
aov4<- Environmental %>%
filter(pH=="ambient")
summary(aov(pHspec ~ as.factor(Header), data=aov4)) #not different between headers
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Header) 3 0.00394 0.0013127 1.906 0.129
## Residuals 305 0.21001 0.0006885
## 251 observations deleted due to missingness
##TANKS
#temp between tanks in high temp
aov1<- Environmental %>%
filter(Temperature=="heated")
summary(aov(Temp ~ as.factor(TankID), data=aov1)) #not different between heated tanks
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(TankID) 11 4.5 0.4077 0.165 0.999
## Residuals 273 675.9 2.4756
## 279 observations deleted due to missingness
#temp between tanks in ambient temp
aov2<- Environmental %>%
filter(Temperature=="ambient")
summary(aov(Temp ~ as.factor(TankID), data=aov2)) #not different between ambient tanks
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(TankID) 11 0.8 0.0687 0.025 1
## Residuals 276 759.8 2.7531
## 276 observations deleted due to missingness
#pH between tanks in acidified
aov3<- Environmental %>%
filter(pH=="acidified")
summary(aov(pHspec ~ as.factor(TankID), data=aov3)) #not different between tanks in acidified
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(TankID) 11 0.0252 0.002288 0.521 0.888
## Residuals 216 0.9495 0.004396
## 340 observations deleted due to missingness
#pH between headers in ambient pH
aov4<- Environmental %>%
filter(pH=="ambient")
summary(aov(pHspec ~ as.factor(TankID), data=aov4)) #not different between headers
## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(TankID) 11 0.00516 0.0004686 0.765 0.675
## Residuals 213 0.13045 0.0006124
## 335 observations deleted due to missingness
tempsummary<-results %>%
select("Day", "Temperature", "Temp") %>%
drop_na(Temp) %>%
group_by(Day, Temperature) %>%
mutate(meanT = mean(Temp)) %>%
group_by(Temperature) %>%
mutate(sd = sd(Temp)) %>%
mutate(se = sd/sqrt(12))
tempplot<-ggplot(data=tempsummary, aes(Day, meanT, color = Temperature)) +
geom_point(size = 2.5, show.legend=FALSE) +
geom_errorbar(aes(ymin = meanT-se, ymax = meanT+se)) +
geom_line(size=1.2) +
scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
scale_y_continuous(name="Temperature (°C)", breaks = seq(20,30, by = 1)) +
scale_shape_discrete(name = NULL,
labels = c("Ambient", "High"))+
scale_color_manual(name = NULL,
values = c("gray40", "firebrick3"),
labels = c("Ambient", "High")) +
ggtitle("")+
theme_classic() +
theme(plot.margin=unit(c(0.3,0.6,0.3,0.3), "cm"))+
theme(legend.margin = margin(0),
legend.position = c(.98,.8),
legend.justification = c("right", "top"),
legend.background = element_blank(),
axis.title = element_text(size = 16, face="bold"),
axis.text = element_text(size = 14),
plot.title = element_text(size=16, face="bold"),
legend.text = element_text(size=14, face="bold"));tempplot
#ggsave(filename="Figures/20200611/TempFig1.png", plot=tempplot, dpi=300, width=7, height=5, units="in")
pHsummary<-results %>%
select("Day","pH", "pHspec") %>%
drop_na(pHspec) %>%
group_by(Day, pH) %>%
mutate(meanpH = mean(pHspec)) %>%
group_by(pH, meanpH) %>%
mutate(sd = sd(pHspec)) %>%
mutate(se = sd/sqrt(12))
pHplot<-ggplot(data=pHsummary, aes(Day, meanpH, color = pH)) +
geom_point(size = 2.5, show.legend=FALSE) +
geom_errorbar(aes(ymin = meanpH-sd, ymax = meanpH+sd)) +
geom_line(size=1.2) +
scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
scale_y_continuous(name="pH (Total Scale)", breaks = seq(7.5,8.1, by = .1)) +
scale_shape_discrete(name = NULL,
labels = c("Low pH", "High pH"))+
scale_color_manual(name = NULL,
values = c("skyblue3", "gray40"),
labels = c("Acidified", "Ambient"),
guide = guide_legend(reverse = TRUE)) +
ggtitle("")+
theme_classic() +
theme(plot.margin=unit(c(0.3,0.6,0.3,0.3), "cm"))+
theme(legend.margin = margin(0),
legend.position = c(.98,.8),
legend.justification = c("right", "top"),
legend.background = element_blank(),
axis.title = element_text(size = 16, face="bold"),
axis.text = element_text(size = 14),
plot.title = element_text(size=16, face="bold"),
legend.text = element_text(size=14, face="bold"));pHplot
#ggsave(filename="Figures/20200611/pHFig1.png", plot=pHplot, dpi=300, width=7, height=5, units="in")
Assemble data for diameters of initial collection of urchins from hatchery
##Initial Collection = Day -24
InitialDiameter <-
#seperate out the initial sizes of urchins on day -24 after initial collection.
results %>%
select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "-24") %>%
#sizes on day -24
select("TankID", "Diameter", "Temperature", "pH", "Treatment")
InitialMean <- mean(InitialDiameter$Diameter)
InitialSE <- std.error(InitialDiameter$Diameter)
InitialMean
## [1] 7.529348
InitialSE
## [1] 0.1068257
Results of linear mixed model effect of diameter (day -24)
InitialMod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=InitialDiameter)
anova(InitialMod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.0005035 0.0005035 1 19 0.1144 0.7389
## pH 0.0102269 0.0102269 1 19 2.3232 0.1439
## Temperature:pH 0.0003840 0.0003840 1 19 0.0872 0.7709
Check assumptions of diameter model (day -24)
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(InitialMod))
## [1] 27 4
shapiro.test(residuals(InitialMod)) #no pass
##
## Shapiro-Wilk normality test
##
## data: residuals(InitialMod)
## W = 0.52108, p-value = 4.878e-11
hist(residuals(InitialMod)) #eh
# 2. Equal variances
leveneTest(residuals(InitialMod)~InitialDiameter$Treatment) #No pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.9592 0.04307 *
## 42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(InitialMod),resid(InitialMod,type="pearson"),col="blue")
Assemble Data for diameters on day 1 of experiment when desired conditions were reached
##After ramp up and desired conditions are reached to begin experiment = Day 1
Day1Diameter <-
#seperate out the initial sizes of urchins on day 1, when desired conditions were met.
results %>%
select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
filter(TankID != "8t" ) %>%
#remove urchin number 8 which died halfway through
drop_na(Diam1, Diam2) %>%
gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>%
#gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
filter(Day == "1") %>%
#sizes on day 1
select("TankID", "Diameter", "Temperature", "pH", "Treatment")
Day1Mean <- mean(Day1Diameter$Diameter)
Day1SE <- std.error(Day1Diameter$Diameter)
Results of linear mixed model effect on diameter (day 1)
Day1Mod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=Day1Diameter)
anova(Day1Mod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.04177 0.04177 1 19 0.1922 0.66606
## pH 0.38094 0.38094 1 19 1.7525 0.20128
## Temperature:pH 0.89969 0.89969 1 19 4.1389 0.05611 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions of diameter model (Day 1)
#check assumptions
# 1. Normality of residuals
qqPlot(residuals(Day1Mod)) #not super normal but... see below
## [1] 1 24
shapiro.test(residuals(Day1Mod)) #passes
##
## Shapiro-Wilk normality test
##
## data: residuals(Day1Mod)
## W = 0.99584, p-value = 0.9998
hist(residuals(Day1Mod)) #looks normalish
# 2. Equal variances
leveneTest(residuals(Day1Mod)~InitialDiameter$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.7309 0.5394
## 42
plot(fitted(Day1Mod),resid(Day1Mod,type="pearson"),col="blue")
Assemble data and calculate means using day 1 as initial diameter and calculate means
## [1] 16.03
## [1] 0.2550688
## # A tibble: 4 x 4
## # Groups: Temperature [2]
## Temperature pH mean s.e.
## <fct> <fct> <dbl> <dbl>
## 1 ambient acidified 327. 16.3
## 2 ambient ambient 323. 11.0
## 3 heated acidified 352. 15.0
## 4 heated ambient 376. 7.15
Growth linear mixed model results
#LMM for growth:
modGrow <-
resultsGrow %>%
lmer(Growth~ Temperature * pH + (1|TankID), data=.)
#Need to decide on model type that we want to use.
summary(modGrow)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ Temperature * pH + (1 | TankID)
## Data: .
##
## REML criterion at convergence: 417.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.90063 -0.32406 0.08141 0.32903 1.94882
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 1929.2 43.92
## Residual 283.1 16.83
## Number of obs: 46, groups: TankID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 326.595 18.578 19.000 17.580 3.28e-13
## Temperatureheated 25.439 26.273 19.000 0.968 0.345
## pHambient -3.424 26.273 19.000 -0.130 0.898
## Temperatureheated:pHambient 27.731 38.073 19.000 0.728 0.475
##
## (Intercept) ***
## Temperatureheated
## pHambient
## Temperatureheated:pHambient
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.707
## pHambient -0.707 0.500
## Tmprtrhtd:H 0.488 -0.690 -0.690
anova(modGrow, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 1169.41 1169.41 1 19 4.1304 0.05634 .
## pH 74.92 74.92 1 19 0.2646 0.61290
## Temperature:pH 150.20 150.20 1 19 0.5305 0.47527
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for growth model
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow))
## [1] 43 20
shapiro.test(residuals(modGrow)) #passes
##
## Shapiro-Wilk normality test
##
## data: residuals(modGrow)
## W = 0.95491, p-value = 0.07263
hist(residuals(modGrow)) #looks normal
# 2. Equal variances
leveneTest(residuals(modGrow)~resultsGrow$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.1906 0.1033
## 42
plot(fitted(modGrow),resid(modGrow,type="pearson"),col="blue")
Assemble data for chosen spine tips (distal end)
### Model 1: calcification at the tips of spines
resultsTip <-
#create data using only the SEM images of the tips of spines.
results %>%
select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>%
#selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
drop_na(RatioSEM) %>%
filter(Chosen == "yes", PartOfSpine == "Tip")
#filters only for tips and chosen side. No dusty images to be analyzed.
Results for spine tips with linear mixed model
#LMM for calcification at the tips of spines
modTip <-
resultsTip %>%
lmer(RatioSEM ~ Temperature * pH + (1|TankID), data = .)
## boundary (singular) fit: see ?isSingular
summary(modTip)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
## Data: .
##
## REML criterion at convergence: 64.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76413 -0.74507 -0.02052 0.63045 2.07630
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 0.0000 0.0000
## Residual 0.1356 0.3682
## Number of obs: 67, groups: TankID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.69118 0.08930 63.00000 18.938 <2e-16
## Temperatureheated -0.06662 0.12453 63.00000 -0.535 0.5945
## pHambient -0.07462 0.12453 63.00000 -0.599 0.5512
## Temperatureheated:pHambient 0.30557 0.18089 63.00000 1.689 0.0961
##
## (Intercept) ***
## Temperatureheated
## pHambient
## Temperatureheated:pHambient .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.717
## pHambient -0.717 0.514
## Tmprtrhtd:H 0.494 -0.688 -0.688
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(modTip, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.10158 0.10158 1 63 0.7492 0.39000
## pH 0.08185 0.08185 1 63 0.6038 0.44006
## Temperature:pH 0.38685 0.38685 1 63 2.8534 0.09612 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for calcification at tip model.
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modTip)) #Normal
## [1] 38 8
shapiro.test(residuals(modTip)) #Pass
##
## Shapiro-Wilk normality test
##
## data: residuals(modTip)
## W = 0.97737, p-value = 0.2605
hist(residuals(modTip))
# 2. Equal variances
leveneTest(residuals(modTip)~resultsTip$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.4545 0.715
## 63
plot(fitted(modTip),resid(modTip,type="pearson"),col="blue")
Conduct post hoc analysis on calcification at tip.
## POST HOC ANALYSIS
emmTip <- emmeans(modTip, ~Temperature*pH, adjust = "tukey") #Estimated marginal means (Least-squares means)
multcomp::cld(emmTip) #create a compact letter display of all pair-wise comparisons
## Temperature pH emmean SE df lower.CL upper.CL .group
## ambient ambient 1.62 0.0868 18.2 1.43 1.80 1
## heated acidified 1.62 0.0868 18.2 1.44 1.81 1
## ambient acidified 1.69 0.0900 18.2 1.50 1.88 1
## heated ambient 1.86 0.0993 16.6 1.65 2.07 1
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
pwpp(emmTip) #Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.
Assemble data for chosen spine bases (proximal end)
### Model 2: calcification in the base of the spines
resultsBase <-
#create data using only the SEM images of the base of spines.
results %>%
select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>%
#selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
drop_na(RatioSEM) %>%
filter(Chosen == "yes", PartOfSpine == "Base")
#filters only for tips and chosen side. No dusty images to be analyzed.
Results for spine bases with linear mixed model
#LMM for calcification at the tips of spines
modBase <-
resultsBase %>%
lmer(RatioSEM~ Temperature * pH + (1|TankID), data=.)
summary(modBase)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
## Data: .
##
## REML criterion at convergence: 98.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5535 -0.6315 -0.1529 0.6918 1.9937
##
## Random effects:
## Groups Name Variance Std.Dev.
## TankID (Intercept) 0.1303 0.3610
## Residual 0.1597 0.3996
## Number of obs: 68, groups: TankID, 23
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.1440 0.1749 18.8077 12.257 2.06e-10
## Temperatureheated 0.3566 0.2487 19.1707 1.434 0.1677
## pHambient 0.6524 0.2474 18.8077 2.637 0.0163
## Temperatureheated:pHambient -0.2224 0.3594 18.9804 -0.619 0.5434
##
## (Intercept) ***
## Temperatureheated
## pHambient *
## Temperatureheated:pHambient
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.703
## pHambient -0.707 0.497
## Tmprtrhtd:H 0.487 -0.692 -0.688
anova(modBase, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature 0.30743 0.30743 1 18.992 1.9251 0.181362
## pH 1.48873 1.48873 1 18.960 9.3223 0.006553 **
## Temperature:pH 0.06114 0.06114 1 18.980 0.3829 0.543424
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for calcification at base model.
## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modBase)) #Normal
## [1] 33 37
shapiro.test(residuals(modBase)) # Pass
##
## Shapiro-Wilk normality test
##
## data: residuals(modBase)
## W = 0.97304, p-value = 0.1472
hist(residuals(modBase))
# 2. Equal variances
leveneTest(residuals(modBase)~resultsBase$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.068 0.369
## 64
plot(fitted(modBase),resid(modBase,type="pearson"),col="blue")
Conduct post hoc analysis on calcification at base
## POST HOC ANALYSIS
emmBase <- emmeans(modBase, ~Temperature*pH, adjust = "tukey")
#Estimated marginal means (Least-squares means)
multcomp::cld(emmBase)
## Temperature pH emmean SE df lower.CL upper.CL .group
## ambient acidified 2.14 0.175 18.8 1.78 2.51 1
## heated acidified 2.50 0.177 19.5 2.13 2.87 12
## ambient ambient 2.80 0.175 18.8 2.43 3.16 12
## heated ambient 2.93 0.192 18.8 2.53 3.33 2
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05
#create a compact letter display of all pair-wise comparisons
pwpp(emmBase)
#Constructs a plot of P values associated with pairwise comparisons of estimated marginal means.
Make Interaction plots for calcification at the base and the tip
Combine both interaction plots at base and tip to be in horizontal line
#Make plots of calcification at tip and base in horizontal line with proper labeling of a and b
#removed legend on the tip, so they can share the same one. But now they are different sizes...
calcificationfig1<-plot_grid(tipplot, baseplot, labels = c("a", "b"), ncol=2, nrow=1, rel_heights= c(1,1), rel_widths = c(0.8,1), label_size = 20, label_y=1, align="h");calcificationfig1
#ggsave(filename="Figures/20200611/CalcificationFig1.png", plot=calcificationfig1, dpi=300, width=16, height=6, units="in")
Assemble data for dropped spines
## Number of loose spines counted at the bottom of tanks at the end of the experiment. These were either broken or shed. Note: these spines were not collected, just counted, so can not determine whether they were broken or shed.
resultsDropped <-
#create data using only the needed variables.
results %>%
select ("Treatment", "Temperature", "pH", "SpineCount", "TankID") %>%
drop_na(SpineCount)
Analyse dropped spines with linear mixed effect model.
##LM of dropped spines - can't make it LMM because error: number of levels of each grouping factor must be < number of observations - only have one count for each TankID.
modDropped <-
resultsDropped %>%
lm(SpineCount~ Temperature * pH, data=.)
summary(modDropped)
##
## Call:
## lm(formula = SpineCount ~ Temperature * pH, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.0000 -6.0833 -0.1667 0.4000 20.0000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.167 3.640 5.814 1.33e-05 ***
## Temperatureheated -1.167 5.148 -0.227 0.82315
## pHambient -21.000 5.148 -4.079 0.00064 ***
## Temperatureheated:pHambient 1.600 7.461 0.214 0.83248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.917 on 19 degrees of freedom
## Multiple R-squared: 0.6088, Adjusted R-squared: 0.547
## F-statistic: 9.855 on 3 and 19 DF, p-value: 0.0003899
anova(modDropped) #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Analysis of Variance Table
##
## Response: SpineCount
## Df Sum Sq Mean Sq F value Pr(>F)
## Temperature 1 1.52 1.52 0.0192 0.8914
## pH 1 2345.78 2345.78 29.4995 3.062e-05 ***
## Temperature:pH 1 3.66 3.66 0.0460 0.8325
## Residuals 19 1510.87 79.52
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Check assumptions for model for dropped spines
##ASSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modDropped)) #not normal
## 1458 1460
## 2 4
shapiro.test(residuals(modDropped)) #faaaiiilll
##
## Shapiro-Wilk normality test
##
## data: residuals(modDropped)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped))
# 2. Equal variances
leveneTest(residuals(modDropped)~resultsDropped$Treatment) # barely pass
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.8487 0.06483 .
## 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modDropped),resid(modDropped,type="pearson"),col="blue")
Analysis for dropped spines does not meet assumptions. Attempt a square root transformation.
##DOES NOT MEET ASSUMPTIONS so ATTEMPT AT Transformation:
resultsDropped$tdata<-(resultsDropped$SpineCount)^1/2
# transformation
modDropped1 <- lm(tdata~ Temperature * pH, data=resultsDropped)
# run lm model of transformed data
###ASSUMPTIONS for Transformation
# 1. Normality of residuals
qqPlot(residuals(modDropped1)) #not normal
## 1458 1460
## 2 4
shapiro.test(residuals(modDropped1)) #fail
##
## Shapiro-Wilk normality test
##
## data: residuals(modDropped1)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped1))
# 2. equal variances
bartlett.test(residuals(modDropped1)~resultsDropped$Treatment) #fail
##
## Bartlett test of homogeneity of variances
##
## data: residuals(modDropped1) by resultsDropped$Treatment
## Bartlett's K-squared = 43.282, df = 3, p-value = 2.144e-09
plot(fitted(modDropped1),resid(modDropped1,type="pearson"),col="blue")
summary(modDropped1)
##
## Call:
## lm(formula = tdata ~ Temperature * pH, data = resultsDropped)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.0000 -3.0417 -0.0833 0.2000 10.0000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.5833 1.8202 5.814 1.33e-05 ***
## Temperatureheated -0.5833 2.5742 -0.227 0.82315
## pHambient -10.5000 2.5742 -4.079 0.00064 ***
## Temperatureheated:pHambient 0.8000 3.7304 0.214 0.83248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.459 on 19 degrees of freedom
## Multiple R-squared: 0.6088, Adjusted R-squared: 0.547
## F-statistic: 9.855 on 3 and 19 DF, p-value: 0.0003899
Anova(modDropped1, type=2) #test for significance of model
## Anova Table (Type II tests)
##
## Response: tdata
## Sum Sq Df F value Pr(>F)
## Temperature 0.23 1 0.0118 0.9146
## pH 586.44 1 29.4995 3.062e-05 ***
## Temperature:pH 0.91 1 0.0460 0.8325
## Residuals 377.72 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Transfomraiton (of any usual kind) of dropped spines was not successful, so we use a non parametric Kruskal Wallis Test.
### Transformations not work, so use non parametric:
kruskal.test(resultsDropped$SpineCount~resultsDropped$Treatment)# SIG*** (p=0.0005563)
##
## Kruskal-Wallis rank sum test
##
## data: resultsDropped$SpineCount by resultsDropped$Treatment
## Kruskal-Wallis chi-squared = 17.505, df = 3, p-value = 0.0005563
Conduct Dunn Post Hoc Test for dropped spines
## POST HOC Pairwise
dunn <-
resultsDropped %>%
dunnTest(SpineCount ~ Treatment,
method = "holm", kw = TRUE, data = .)
dunnph <- dunn$res
cldList(P.adj ~ Comparison, data = dunnph, threshold = 0.05)
## Group Letter MonoLetter
## 1 AMB a a
## 2 OA b b
## 3 T ac a c
## 4 T/OA bc bc
droppedsummary<-results %>%
select("Temperature","pH","SpineCount") %>%
drop_na(SpineCount) %>%
group_by(pH, Temperature) %>%
summarize(mean = mean(SpineCount), se = std.error(SpineCount), N = length(SpineCount))
dropplot<-ggplot(data=droppedsummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
geom_point(size = 3, position=position_dodge(0.1)) +
geom_line(aes(group=pH), position=position_dodge(0.1)) +
xlab("Temperature Treatment") +
ylab("Dropped Spines")+
geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
ylim(0,30) +
theme_bw() +
scale_x_discrete(limits=c("ambient", "heated"),
label=c("Ambient", "High"))+
scale_colour_manual(name = "pH Treatment",
values = c("skyblue3", "gray40"),
labels = c("Acidified", "Ambient"),
guide = guide_legend(reverse = TRUE)) +
#geom_text(x=1.5, y=11, label="p(pH)<0.001", size=6, color="black") +
#geom_text(x=1.5, y=9, label="p(Temperature)=0.91", size=6, color="darkgray") +
#geom_text(x=1.5, y=7, label="p(Interaction)=0.83", size=6, color="darkgray") +
theme_classic()+
theme(legend.position = "right",
legend.background = element_blank(),
axis.title = element_text(size = 16, face="bold"),
axis.text = element_text(size = 14),
plot.title = element_text(size=16, face="bold"),
legend.text = element_text(size=14),
legend.title = element_text(size=16, face="bold"));dropplot
## Warning: Removed 1 rows containing missing values (geom_errorbar).
#ggsave(filename="Figures/20200611/DroppedSpinesFig.png", plot=dropplot, dpi=300, width=8, height=6, units="in")